function S = ECI2CR3BP(S)
u = [398600.44; 4902.80007; 0];
    %[km^3/s^2]Gravitational parameters of the Earth, Moon, and satellite.
    
    r12 = 384400;
    %[km]Distance between the Earth and the Moon.
    
    rho(1) = u(1) / sum(u);
    %[]Ratio of the gravitational parameters.
    
    S(1:9,:) = S(1:9,:);
    %[]Normalizes all state vectors.
    
    Wcm = sqrt(sum(u)/r12^3)*[0;0;1];
    %[] angular velocity in Cr3bp frame
    
    n = size(S,2);
    %[]Number of position vectors.
    
    Rm = S(4:6,:);
    %[]Moon positions WRT the CM in ECI coordinates.
    Rcme = -S(1:3,:);
    Vm = S(13:15,:);
    %[1/s]Moon velocities WRT the CM in ECI coordinates.
    
    Rsat = S(7:9,:);
    %[]Satellite positions WRT the CM in ECI coordinates.
    Vsat = S(16:18,:);
    % Satellite velocity wrt CM in ECI
    
    for k = 1:n
        
        Hm = cross(Rm(:,k),Vm(:,k));
        %[1/s]Current specific angular momentum of the Moon WRT the CM in ECI coordinates.
        
        Rhat = Rm(:,k) / norm(Rm(:,k));
        %[]Current CR3BP radial direction in ECI coordinates.
        
        Nhat = Hm / norm(Hm);
        %[]Current CR3BP normal direction in ECI coordinates.
        
        That = cross(Nhat,Rhat);
        %[]Current CR3BP tangential direction in ECI coordinates.
        
        Cr3bpToEci = [Rhat, That, Nhat];
        %[]Matrix that transforms vectors from CR3BP to ECI coordinates.
        
        Rsat(:,k) = transpose(Cr3bpToEci) * Rsat(:,k);
        %[]Current satellite position WRT the CM in ECI coordinates.
        
        Rcme(:,k) = transpose(Cr3bpToEci) * Rcme(:,k);
        
        Vsat(:,k) = transpose(Cr3bpToEci) * Vsat(:,k)-cross(Wcm,Rcme(:,k)) - cross(Wcm,Rsat(:,k));
    end

S = [Rsat;Vsat];
end